library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(bayesplot)
## This is bayesplot version 1.7.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
library(cowplot)
## 
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
##   default ggplot2 theme anymore. To recover the previous
##   behavior, execute:
##   theme_set(theme_cowplot())
## ********************************************************
library(testthat)
## 
## Attaching package: 'testthat'
## The following object is masked from 'package:dplyr':
## 
##     matches
library(parallel)
library(purrr)
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:testthat':
## 
##     is_null
rstan_options(auto_write=TRUE)
# generate samples from a truncated normal distribution. n = how many samples, mean and sd are mean and sd, and min and max are the limits of the distribution/truncation points. possibly my first ever use of while loops. Don't @ me.
rtnorm <- function(n, mean, sd, min, max) {
    x <- rnorm(n, mean=mean, sd=sd)
    x <- x[x >= min & x <= max]
    while(length(x) < n) {
        newx <- rnorm(1, mean=mean, sd=sd)
        while(newx <= min | newx >= max) {
            newx <- rnorm(1, mean=mean, sd=sd)
        }
        x <- c(x, newx)
    }
    length(x)==n
    return(x)
}

# simulate data from an ordinal logistic model and format it as input for stan. input is a list for stan, as is output.
simulate_data <- function(input) {
    simu <- rstan::stan(file='dirichlet prior/covar_sim.stan', iter=1, chains=1, algorithm="Fixed_param", data=input)

    simu_params <- rstan::extract(simu)

    input_data_for_model <- list("N" = input$N, "K" = input$K, "x" = input$x, "y" = array(simu_params$y[1,]))
    return(input_data_for_model)
}

# plot simulated data for a sanity check.
simplot <- function(datalist) {
    inputdf <- data.frame(datalist)
    p1 <- ggplot(inputdf, aes(x=x, y=y)) +
        geom_jitter(shape=1, height=0.2) +
        geom_vline(xintercept = h) +
        ggtitle("Simulated data with cutpoints")
    print(p1)
    p2 <- ggplot(inputdf, aes(x=x, colour=as.factor(y))) +
        stat_ecdf() +
        geom_vline(xintercept=h) +
        theme(legend.position = "none") +
        ggtitle("Cumulative x for 3 states")
    print(p2)
    p3 <- cowplot::plot_grid(p1,p2)
    print(p3)
}

# append the correct shape (gamma prior on cutpoints) and rate (exponential prior on beta) parameters to the simulated data list and fit a model with a gamma prior. simdat is a list of datasets simulated from stan models and pars is a vector of additional parameters. Relies on global params beta_slow and beta_fast.
fit_gamma_model <- function(simdatlist, pars) {
    #choose whether to use data simulated with a rapid or slow transition
    if (pars$beta == beta_slow) {
        simdat <- simdatlist$simdat_slow
    }
    if (pars$beta == beta_fast) {
        simdat <- simdatlist$simdat_fast
    }
    #extract parameters for prior distribtuions
    simdat$shape <- pars$shape
    simdat$rate <- pars$beta_rate

    #fit the model
    fitgam <- stan(file='dirichlet prior/gamma/gamma_covar.stan', data=simdat, chains=4)
    return(fitgam)
}

## append a label (string) to all columnnames in a dataframe (x)
label_names <- function(x, label) {
    colnames(x) <- paste0(colnames(x), "_", label)
    return(x)
}

# function to plot modeled parameters with label being a string to label the graph (usually prior and whether groups were included and pardf is the modeled parameter dataframe) trueparams is a one row dataframe of true parameters. h is a global parameter have fun!
parplot <- function(pars) {
    cutplot <- mcmc_areas(pars, regex_pars="c.\\d_model") +
        geom_vline(xintercept=c(pars$c.1_true, pars$c.2_true))

    betaplot <- mcmc_areas(pars, pars="beta_model") +
        geom_vline(xintercept=pars$beta_true)

    h1plot <- mcmc_areas(pars, regex_pars = "h1") +
        geom_vline(xintercept=h[1])
    h2plot <- mcmc_areas(pars, regex_pars = "h2") +
        geom_vline(xintercept=h[2])
    allplot <- cowplot::plot_grid(cutplot, betaplot, h1plot, h2plot,
                                  nrow=1, ncol=4,
                                  labels=paste("model", unique(pars$modelid_true),
                                               "rate=", unique(pars$beta_rate_true),
                                               "shape=", unique(pars$shape_true)))
    print(allplot)
}

# function to calculate the difference between modeled parameters (in the model_params dataframe) and true parameters (h is globally declared). Where true params is a dataframe
posterior_differencer <- function(pars) {
    c1_diff <- pars$c.1_model - pars$c.1_true
    c2_diff <- pars$c.2_model - pars$c.1_true
    h1_diff <- pars$h1_model - h[1]
    h2_diff <- pars$h2_model - h[2]
    beta_diff <- pars$beta_model - pars$beta_true
    diffframe <- data.frame(c1_diff, c2_diff, h1_diff, h2_diff, beta_diff)
    return(diffframe)
}

# function to plot histograms of differences between true params and modeled params (model_params dataframe).
diffplotter <- function(diffs, pars) {
    #diffs <- posterior_differencer(pars)
    cuts <- mcmc_areas(diffs, regex_pars = "c") +
        ggtitle("", subtitle = "differences between modeled and true params")+
        xlim(c(-15,15))
    opars <- mcmc_areas(diffs, pars=c("h1_diff", "h2_diff", "beta_diff")) +
        xlim(c(-1,1))
    cowplot::plot_grid(cuts, opars, labels=paste("model", unique(pars$modelid_true),
                                                 "rate=", unique(pars$beta_rate_true),
                                                 "shape=", unique(pars$shape_true)))
}

Prior choice

In my original phenology model, I use

Michael Betancourt thinks that an induced dirichlet prior may be a good choice for ordinal logistic models.

Normally, priors are bottom up - you choose a distribution for each prior. This is tricky for the cutpoints in an ordinal logistic model because they’re defined on an abstract latent space so you can’t easily use your domain expertise, but a good prior is incredibly important in these models, because ordered logistic models with covariates are inherently non-identifiable. (Because \(beta\)s and cutpoints depend on one another.)

Betancourt says that > To avoid the non-identifiability of the interior cut points from propagating to the posterior distribution we need a principled prior model that can exploit domain expertise to consistently regularize all of the internal cut points at the same time. Regularization of the cut points, however, is subtle given their ordering constraint. Moreover domain expertise is awkward to apply on the abstract latent space where the internal cut points are defined.

Goals

  1. Understand the dynamics of the gamma and induced dirichlet priors and 2) compare them

Simulate data

First let’s simulate some data that’s kind of like mine.

We have three potential states K and a latent effect/covariate x. Data is most likely to be collected around state 2. x is always positive.

I’ll simulate 2 datasets - one with a fast transition speed (\(\beta = 2\)) and one with a slow transition speed (\(\beta = 0.5\)). I want the transitions to occur at the same x for both datasets. The halfway transition points h will be at x= 8 and x=12. So the cutpoints c for the slow transition will be at 4 and 6 and for the fast transition at 16 and 24.

Simulate datasets

N <- 500
K <- 3

beta_slow <- 0.5
beta_fast <- 2

cutpoints_slow <- c(4,6)
cutpoints_fast <- c(16,24)

h <- cutpoints_slow/beta_slow #half transition points, engineered to be identical for slow and fast transitions
x <- rtnorm(n=N, mean=mean(h), sd=2, min=0, max=20) #covariate

testthat::test_that("half transition points identical", {
  testthat::expect_equal(cutpoints_slow/beta_slow, cutpoints_fast/beta_fast)
})
inputs_for_sim_slow <- list("N" = N, "K" = K, "c" = cutpoints_slow, "beta"=beta_slow, "h" = h, "x" = x)
inputs_for_sim_fast <- list("N" = N, "K" = K, "c" = cutpoints_fast, "beta"=beta_fast, "h" = h, "x" = x)

# simulate data, graph it, and prepare data for model fitting

simdat_slow <- simulate_data(input=inputs_for_sim_slow)
## 
## SAMPLING FOR MODEL 'covar_sim' NOW (CHAIN 1).
## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
## Chain 1:                0.000153 seconds (Sampling)
## Chain 1:                0.000153 seconds (Total)
## Chain 1:
slowsimplot <- simplot(simdat_slow)

simdat_fast <- simulate_data(input=inputs_for_sim_fast)
## 
## SAMPLING FOR MODEL 'covar_sim' NOW (CHAIN 1).
## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
## Chain 1:                0.00015 seconds (Sampling)
## Chain 1:                0.00015 seconds (Total)
## Chain 1:
fastsimplot <- simplot(simdat_fast)

cowplot::plot_grid(slowsimplot, fastsimplot, nrow=2, labels=c("slow", "fast"), vjust=2)

Recapture parameters: Gamma

Now I want to try to recapture parameters.

I need to know how good I have to be at choosing the shape parameter for the gamma prior on the cutpoints to recapture parameters as well. So I’ll choose a prior for gamma that roughly centers it between the cutpoints, which I think is a “good” option (though given the long right tail, maybe a lower value would be better). Then halve and double it.

I also want to understand how the exponential prior on beta affects the ability to recapture parameters. I’ll try rates of 1,2,3.

Gamma prior with covariate

beta_rate <- c(1:3) # rate parameters for exponential prior on beta

# shape parameters for gamma prior on cutpoints, for slow and fast transitions
shape_slow <- mean(cutpoints_slow)
shape_slow_small <- shape_slow/2
shape_slow_big <- shape_slow*2

slowshapes <- c(shape_slow, shape_slow_small, shape_slow_big)

shape_fast <- mean(cutpoints_fast)
shape_fast_small <- shape_fast/2
shape_fast_big <- shape_fast*2

fastshapes <- c(shape_fast, shape_fast_big, shape_fast_small)

slowframe <- data.frame(beta = beta_slow, shape = slowshapes, beta_rate=beta_rate, c.1 = cutpoints_slow[1], c.2=cutpoints_slow[2], h1=h[1], h2=h[2])
fastframe <- data.frame(beta=beta_fast, shape= fastshapes, beta_rate=beta_rate, c.1 = cutpoints_fast[1], c.2=cutpoints_fast[2], h1=h[1], h2=h[2])
parframe <- rbind(slowframe, fastframe) %>%
  tidyr::expand(tidyr::nesting(beta, c.1, c.2, shape), h1, h2, beta_rate)
parframe$modelid <- 1:nrow(parframe) #label the models
parframe <- dplyr::select(parframe, modelid, beta, c.1, c.2, h1, h2, shape, beta_rate)

knitr::kable(parframe, caption="beta parameters used to simulate data for a slow (0.5) and fast (2) transion. And shape and rate parameters to be used in the priors on beta and cutpoints to attempt to recover parameters")
beta parameters used to simulate data for a slow (0.5) and fast (2) transion. And shape and rate parameters to be used in the priors on beta and cutpoints to attempt to recover parameters
modelid beta c.1 c.2 h1 h2 shape beta_rate
1 0.5 4 6 8 12 2.5 1
2 0.5 4 6 8 12 2.5 2
3 0.5 4 6 8 12 2.5 3
4 0.5 4 6 8 12 5.0 1
5 0.5 4 6 8 12 5.0 2
6 0.5 4 6 8 12 5.0 3
7 0.5 4 6 8 12 10.0 1
8 0.5 4 6 8 12 10.0 2
9 0.5 4 6 8 12 10.0 3
10 2.0 16 24 8 12 10.0 1
11 2.0 16 24 8 12 10.0 2
12 2.0 16 24 8 12 10.0 3
13 2.0 16 24 8 12 20.0 1
14 2.0 16 24 8 12 20.0 2
15 2.0 16 24 8 12 20.0 3
16 2.0 16 24 8 12 40.0 1
17 2.0 16 24 8 12 40.0 2
18 2.0 16 24 8 12 40.0 3

So there are 2 simulated datasets (beta=0.5 or beta=2) and I’m going to try to fit them both with 3 rates for the beta’s exponential prior and 3 shapes for the cutpoints’ gamma prior, a total of 18 model runs.

parlist <- split(parframe, seq(nrow(parframe))) # format parframe so it works with parLapply better
names(parlist) <- paste0("m", 1:nrow(parframe))
datlist <- list(simdat_slow=simdat_slow, simdat_fast=simdat_fast) # list of simulated datasets

# run all models, parallelized

# make a cluster, leaving 20 cores free for other folks
no_cores <- parallel::detectCores() - 20
cl <- parallel::makeCluster(no_cores)

# export the stuff you need to run on the cluster
parallel::clusterExport(cl, c("fit_gamma_model", "parlist", "beta_fast", "beta_slow", "datlist"))
parallel::clusterEvalQ(cl, c(library(rstan), library(StanHeaders)))
## [[1]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[2]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[3]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[4]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[5]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[6]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[7]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[8]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[9]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[10]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[11]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[12]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[13]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[14]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[15]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[16]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[17]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[18]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[19]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## 
## [[20]]
##  [1] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
##  [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [11] "rstan"       "ggplot2"     "StanHeaders" "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"
# fit the models
gammafits <- parallel::parLapply(cl, parlist, function(x) {fit_gamma_model(simdatlist = datlist, pars=x)})

parallel::stopCluster(cl) #close the cluster

# extract params
paramsgam <- lapply(gammafits,
                    function(x) {data.frame(rstan::extract(x) ) } )
#param summary
sum <- lapply(paramsgam, summary)

#bind true and model pars even tho it will make a giant df


paramsgam <- map(paramsgam, label_names, label="model")
parlist <- map(parlist, label_names, label="true")

paramsgam <- map2(paramsgam, parlist, cbind)

#plot params and diffs
map(paramsgam, parplot)

## $m1

## 
## $m2

## 
## $m3

## 
## $m4

## 
## $m5

## 
## $m6

## 
## $m7

## 
## $m8

## 
## $m9

## 
## $m10

## 
## $m11

## 
## $m12

## 
## $m13

## 
## $m14

## 
## $m15

## 
## $m16

## 
## $m17

## 
## $m18

map(paramsgam, posterior_differencer) %>%
    map2(.y=paramsgam, .f=diffplotter)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 2 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_segment).
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## $m1

## 
## $m2

## 
## $m3

## 
## $m4

## 
## $m5

## 
## $m6

## 
## $m7

## 
## $m8

## 
## $m9

## 
## $m10

## 
## $m11

## 
## $m12

## 
## $m13

## 
## $m14

## 
## $m15

## 
## $m16

## 
## $m17

## 
## $m18

# calculate whether true value is in 50% hpdi

# kludge HPDI
HPDIlow <- function(x, prob) {
    HPDI <- rethinking::HPDI(x, prob=prob)
    return(HPDI[1])
}

HPDIhigh <- function(x, prob) {
    HPDI <- rethinking::HPDI(x, prob=prob)
    return(HPDI[2])
}


calc_HPDI <- function(params, prob) {
    low <- params %>% dplyr::summarise_at(vars(ends_with("model")), HPDIlow, prob=0.5)
    high <- params %>% dplyr::summarise_at(vars(ends_with("model")), HPDIhigh, prob=0.5)

    # awkward formatting
    hdpis <- dplyr::full_join(low, high) %>%
        select(-contains("lp"))
    colnames(hdpis) <- stringr::str_replace(colnames(hdpis), "_model", "")

    # true param
    true <- params %>% dplyr::summarise_at(vars(ends_with("true")), unique)
    colnames(true) <- stringr::str_replace(colnames(true), "_true", "")
    true <- select(true, colnames(hdpis))

    # more awkward formatting
    compframe <- dplyr::full_join(hdpis, true) %>%
        t(.) %>%
        data.frame()
    colnames(compframe) <- c("low", "high", "true")
    compframe$params <- rownames(compframe)

    # true param in interval?
    tf <- compframe %>% mutate(inint = true > low & true < high)
    return(tf)

}

in50 <- map(paramsgam, calc_HPDI, prob=0.5)
## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")
## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")
print(in50)
## $m1
##          low       high true params inint
## 1  3.2501236  3.9272662  4.0    c.1 FALSE
## 2  4.9079239  5.6276046  6.0    c.2 FALSE
## 3  0.4066015  0.4748422  0.5   beta FALSE
## 4  7.8773879  8.2708234  8.0     h1  TRUE
## 5 11.7804592 12.1964284 12.0     h2  TRUE
## 
## $m2
##          low       high true params inint
## 1  3.1747532  3.8290328  4.0    c.1 FALSE
## 2  4.9277306  5.6227433  6.0    c.2 FALSE
## 3  0.4031203  0.4692065  0.5   beta FALSE
## 4  7.8930745  8.2767062  8.0     h1  TRUE
## 5 11.8244379 12.2359328 12.0     h2  TRUE
## 
## $m3
##          low       high true params inint
## 1  3.1904716  3.8115469  4.0    c.1 FALSE
## 2  4.9024314  5.5857277  6.0    c.2 FALSE
## 3  0.4086681  0.4723668  0.5   beta FALSE
## 4  7.8722504  8.2621884  8.0     h1  TRUE
## 5 11.8230640 12.2186223 12.0     h2  TRUE
## 
## $m4
##          low       high true params inint
## 1  3.4753182  4.1142259  4.0    c.1  TRUE
## 2  5.2798238  5.9367518  6.0    c.2 FALSE
## 3  0.4363857  0.4991233  0.5   beta FALSE
## 4  8.0089769  8.3662952  8.0     h1 FALSE
## 5 11.7793695 12.1459087 12.0     h2  TRUE
## 
## $m5
##          low       high true params inint
## 1  3.4433367  4.0945892  4.0    c.1  TRUE
## 2  5.1325837  5.8360014  6.0    c.2 FALSE
## 3  0.4246286  0.4904984  0.5   beta FALSE
## 4  8.0189219  8.3836955  8.0     h1 FALSE
## 5 11.7668734 12.1422131 12.0     h2  TRUE
## 
## $m6
##         low       high true params inint
## 1  3.370690  4.0211574  4.0    c.1  TRUE
## 2  5.186302  5.8966790  6.0    c.2 FALSE
## 3  0.418465  0.4850303  0.5   beta FALSE
## 4  7.999701  8.3563466  8.0     h1  TRUE
## 5 11.791000 12.1658031 12.0     h2  TRUE
## 
## $m7
##          low       high true params inint
## 1  3.9334486  4.5523476  4.0    c.1  TRUE
## 2  5.7414360  6.4048148  6.0    c.2  TRUE
## 3  0.4744432  0.5376042  0.5   beta  TRUE
## 4  8.2128315  8.5227315  8.0     h1 FALSE
## 5 11.6305653 11.9513433 12.0     h2 FALSE
## 
## $m8
##          low      high true params inint
## 1  3.9378608  4.557596  4.0    c.1  TRUE
## 2  5.7729105  6.447308  6.0    c.2  TRUE
## 3  0.4702384  0.533441  0.5   beta  TRUE
## 4  8.2315496  8.543833  8.0     h1 FALSE
## 5 11.6784857 11.993660 12.0     h2 FALSE
## 
## $m9
##          low       high true params inint
## 1  3.9731254  4.5753073  4.0    c.1  TRUE
## 2  5.6989545  6.3508008  6.0    c.2  TRUE
## 3  0.4805663  0.5419001  0.5   beta  TRUE
## 4  8.1834663  8.5009195  8.0     h1 FALSE
## 5 11.6774571 11.9984272 12.0     h2 FALSE
## 
## $m10
##         low      high true params inint
## 1 11.159384 12.494009   16    c.1 FALSE
## 2 18.014987 19.826380   24    c.2 FALSE
## 3  1.448803  1.602843    2   beta FALSE
## 4  7.672810  7.841917    8     h1 FALSE
## 5 12.192592 12.356752   12     h2 FALSE
## 
## $m11
##         low      high true params inint
## 1 11.185234 12.402412   16    c.1 FALSE
## 2 17.852498 19.520479   24    c.2 FALSE
## 3  1.445014  1.587621    2   beta FALSE
## 4  7.656411  7.817452    8     h1 FALSE
## 5 12.170768 12.333943   12     h2 FALSE
## 
## $m12
##         low      high true params inint
## 1 10.848001 12.049589   16    c.1 FALSE
## 2 17.692949 19.370861   24    c.2 FALSE
## 3  1.421928  1.564056    2   beta FALSE
## 4  7.664332  7.829884    8     h1 FALSE
## 5 12.197229 12.360644   12     h2 FALSE
## 
## $m13
##         low      high true params inint
## 1 12.524896 13.899434   16    c.1 FALSE
## 2 19.839749 21.731258   24    c.2 FALSE
## 3  1.604186  1.765372    2   beta FALSE
## 4  7.760157  7.913378    8     h1 FALSE
## 5 12.146768 12.297673   12     h2 FALSE
## 
## $m14
##         low      high true params inint
## 1 12.320277 13.634744   16    c.1 FALSE
## 2 19.447600 21.265906   24    c.2 FALSE
## 3  1.592442  1.745851    2   beta FALSE
## 4  7.735040  7.883800    8     h1 FALSE
## 5 12.148636 12.299303   12     h2 FALSE
## 
## $m15
##         low      high true params inint
## 1 12.213129 13.531782   16    c.1 FALSE
## 2 19.267255 21.137336   24    c.2 FALSE
## 3  1.575907  1.733495    2   beta FALSE
## 4  7.729892  7.877785    8     h1 FALSE
## 5 12.143976 12.292382   12     h2 FALSE
## 
## $m16
##         low      high true params inint
## 1 15.299240 16.740585   16    c.1  TRUE
## 2 23.408511 25.483513   24    c.2  TRUE
## 3  1.930559  2.106695    2   beta  TRUE
## 4  7.854474  7.987955    8     h1 FALSE
## 5 12.089570 12.224809   12     h2 FALSE
## 
## $m17
##         low      high true params inint
## 1 15.322033 16.795263   16    c.1  TRUE
## 2 23.497265 25.605068   24    c.2  TRUE
## 3  1.939864  2.115915    2   beta  TRUE
## 4  7.845753  7.980946    8     h1 FALSE
## 5 12.070039 12.206998   12     h2 FALSE
## 
## $m18
##         low      high true params inint
## 1 14.975978 16.394557   16    c.1  TRUE
## 2 23.018115 24.980390   24    c.2  TRUE
## 3  1.907683  2.074526    2   beta  TRUE
## 4  7.856500  7.991180    8     h1 FALSE
## 5 12.098570 12.233131   12     h2 FALSE
in75 <- map(paramsgam, calc_HPDI, prob=0.75)
## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")
print(in75)
## $m1
##          low       high true params inint
## 1  3.2501236  3.9272662  4.0    c.1 FALSE
## 2  4.9079239  5.6276046  6.0    c.2 FALSE
## 3  0.4066015  0.4748422  0.5   beta FALSE
## 4  7.8773879  8.2708234  8.0     h1  TRUE
## 5 11.7804592 12.1964284 12.0     h2  TRUE
## 
## $m2
##          low       high true params inint
## 1  3.1747532  3.8290328  4.0    c.1 FALSE
## 2  4.9277306  5.6227433  6.0    c.2 FALSE
## 3  0.4031203  0.4692065  0.5   beta FALSE
## 4  7.8930745  8.2767062  8.0     h1  TRUE
## 5 11.8244379 12.2359328 12.0     h2  TRUE
## 
## $m3
##          low       high true params inint
## 1  3.1904716  3.8115469  4.0    c.1 FALSE
## 2  4.9024314  5.5857277  6.0    c.2 FALSE
## 3  0.4086681  0.4723668  0.5   beta FALSE
## 4  7.8722504  8.2621884  8.0     h1  TRUE
## 5 11.8230640 12.2186223 12.0     h2  TRUE
## 
## $m4
##          low       high true params inint
## 1  3.4753182  4.1142259  4.0    c.1  TRUE
## 2  5.2798238  5.9367518  6.0    c.2 FALSE
## 3  0.4363857  0.4991233  0.5   beta FALSE
## 4  8.0089769  8.3662952  8.0     h1 FALSE
## 5 11.7793695 12.1459087 12.0     h2  TRUE
## 
## $m5
##          low       high true params inint
## 1  3.4433367  4.0945892  4.0    c.1  TRUE
## 2  5.1325837  5.8360014  6.0    c.2 FALSE
## 3  0.4246286  0.4904984  0.5   beta FALSE
## 4  8.0189219  8.3836955  8.0     h1 FALSE
## 5 11.7668734 12.1422131 12.0     h2  TRUE
## 
## $m6
##         low       high true params inint
## 1  3.370690  4.0211574  4.0    c.1  TRUE
## 2  5.186302  5.8966790  6.0    c.2 FALSE
## 3  0.418465  0.4850303  0.5   beta FALSE
## 4  7.999701  8.3563466  8.0     h1  TRUE
## 5 11.791000 12.1658031 12.0     h2  TRUE
## 
## $m7
##          low       high true params inint
## 1  3.9334486  4.5523476  4.0    c.1  TRUE
## 2  5.7414360  6.4048148  6.0    c.2  TRUE
## 3  0.4744432  0.5376042  0.5   beta  TRUE
## 4  8.2128315  8.5227315  8.0     h1 FALSE
## 5 11.6305653 11.9513433 12.0     h2 FALSE
## 
## $m8
##          low      high true params inint
## 1  3.9378608  4.557596  4.0    c.1  TRUE
## 2  5.7729105  6.447308  6.0    c.2  TRUE
## 3  0.4702384  0.533441  0.5   beta  TRUE
## 4  8.2315496  8.543833  8.0     h1 FALSE
## 5 11.6784857 11.993660 12.0     h2 FALSE
## 
## $m9
##          low       high true params inint
## 1  3.9731254  4.5753073  4.0    c.1  TRUE
## 2  5.6989545  6.3508008  6.0    c.2  TRUE
## 3  0.4805663  0.5419001  0.5   beta  TRUE
## 4  8.1834663  8.5009195  8.0     h1 FALSE
## 5 11.6774571 11.9984272 12.0     h2 FALSE
## 
## $m10
##         low      high true params inint
## 1 11.159384 12.494009   16    c.1 FALSE
## 2 18.014987 19.826380   24    c.2 FALSE
## 3  1.448803  1.602843    2   beta FALSE
## 4  7.672810  7.841917    8     h1 FALSE
## 5 12.192592 12.356752   12     h2 FALSE
## 
## $m11
##         low      high true params inint
## 1 11.185234 12.402412   16    c.1 FALSE
## 2 17.852498 19.520479   24    c.2 FALSE
## 3  1.445014  1.587621    2   beta FALSE
## 4  7.656411  7.817452    8     h1 FALSE
## 5 12.170768 12.333943   12     h2 FALSE
## 
## $m12
##         low      high true params inint
## 1 10.848001 12.049589   16    c.1 FALSE
## 2 17.692949 19.370861   24    c.2 FALSE
## 3  1.421928  1.564056    2   beta FALSE
## 4  7.664332  7.829884    8     h1 FALSE
## 5 12.197229 12.360644   12     h2 FALSE
## 
## $m13
##         low      high true params inint
## 1 12.524896 13.899434   16    c.1 FALSE
## 2 19.839749 21.731258   24    c.2 FALSE
## 3  1.604186  1.765372    2   beta FALSE
## 4  7.760157  7.913378    8     h1 FALSE
## 5 12.146768 12.297673   12     h2 FALSE
## 
## $m14
##         low      high true params inint
## 1 12.320277 13.634744   16    c.1 FALSE
## 2 19.447600 21.265906   24    c.2 FALSE
## 3  1.592442  1.745851    2   beta FALSE
## 4  7.735040  7.883800    8     h1 FALSE
## 5 12.148636 12.299303   12     h2 FALSE
## 
## $m15
##         low      high true params inint
## 1 12.213129 13.531782   16    c.1 FALSE
## 2 19.267255 21.137336   24    c.2 FALSE
## 3  1.575907  1.733495    2   beta FALSE
## 4  7.729892  7.877785    8     h1 FALSE
## 5 12.143976 12.292382   12     h2 FALSE
## 
## $m16
##         low      high true params inint
## 1 15.299240 16.740585   16    c.1  TRUE
## 2 23.408511 25.483513   24    c.2  TRUE
## 3  1.930559  2.106695    2   beta  TRUE
## 4  7.854474  7.987955    8     h1 FALSE
## 5 12.089570 12.224809   12     h2 FALSE
## 
## $m17
##         low      high true params inint
## 1 15.322033 16.795263   16    c.1  TRUE
## 2 23.497265 25.605068   24    c.2  TRUE
## 3  1.939864  2.115915    2   beta  TRUE
## 4  7.845753  7.980946    8     h1 FALSE
## 5 12.070039 12.206998   12     h2 FALSE
## 
## $m18
##         low      high true params inint
## 1 14.975978 16.394557   16    c.1  TRUE
## 2 23.018115 24.980390   24    c.2  TRUE
## 3  1.907683  2.074526    2   beta  TRUE
## 4  7.856500  7.991180    8     h1 FALSE
## 5 12.098570 12.233131   12     h2 FALSE
in90 <- map(paramsgam, calc_HPDI, prob=0.90)
## Joining, by = c("c.1_model", "c.2_model", "beta_model", "h1_model", "h2_model",
## "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")Joining, by = c("c.1_model",
## "c.2_model", "beta_model", "h1_model", "h2_model", "lp___model")
## Warning: Column `c.1_model` has different attributes on LHS and RHS of join
## Warning: Column `c.2_model` has different attributes on LHS and RHS of join
## Warning: Column `beta_model` has different attributes on LHS and RHS of join
## Warning: Column `h1_model` has different attributes on LHS and RHS of join
## Warning: Column `h2_model` has different attributes on LHS and RHS of join
## Warning: Column `lp___model` has different attributes on LHS and RHS of join
## Joining, by = c("c.1", "c.2", "beta", "h1", "h2")
print(in90)
## $m1
##          low       high true params inint
## 1  3.2501236  3.9272662  4.0    c.1 FALSE
## 2  4.9079239  5.6276046  6.0    c.2 FALSE
## 3  0.4066015  0.4748422  0.5   beta FALSE
## 4  7.8773879  8.2708234  8.0     h1  TRUE
## 5 11.7804592 12.1964284 12.0     h2  TRUE
## 
## $m2
##          low       high true params inint
## 1  3.1747532  3.8290328  4.0    c.1 FALSE
## 2  4.9277306  5.6227433  6.0    c.2 FALSE
## 3  0.4031203  0.4692065  0.5   beta FALSE
## 4  7.8930745  8.2767062  8.0     h1  TRUE
## 5 11.8244379 12.2359328 12.0     h2  TRUE
## 
## $m3
##          low       high true params inint
## 1  3.1904716  3.8115469  4.0    c.1 FALSE
## 2  4.9024314  5.5857277  6.0    c.2 FALSE
## 3  0.4086681  0.4723668  0.5   beta FALSE
## 4  7.8722504  8.2621884  8.0     h1  TRUE
## 5 11.8230640 12.2186223 12.0     h2  TRUE
## 
## $m4
##          low       high true params inint
## 1  3.4753182  4.1142259  4.0    c.1  TRUE
## 2  5.2798238  5.9367518  6.0    c.2 FALSE
## 3  0.4363857  0.4991233  0.5   beta FALSE
## 4  8.0089769  8.3662952  8.0     h1 FALSE
## 5 11.7793695 12.1459087 12.0     h2  TRUE
## 
## $m5
##          low       high true params inint
## 1  3.4433367  4.0945892  4.0    c.1  TRUE
## 2  5.1325837  5.8360014  6.0    c.2 FALSE
## 3  0.4246286  0.4904984  0.5   beta FALSE
## 4  8.0189219  8.3836955  8.0     h1 FALSE
## 5 11.7668734 12.1422131 12.0     h2  TRUE
## 
## $m6
##         low       high true params inint
## 1  3.370690  4.0211574  4.0    c.1  TRUE
## 2  5.186302  5.8966790  6.0    c.2 FALSE
## 3  0.418465  0.4850303  0.5   beta FALSE
## 4  7.999701  8.3563466  8.0     h1  TRUE
## 5 11.791000 12.1658031 12.0     h2  TRUE
## 
## $m7
##          low       high true params inint
## 1  3.9334486  4.5523476  4.0    c.1  TRUE
## 2  5.7414360  6.4048148  6.0    c.2  TRUE
## 3  0.4744432  0.5376042  0.5   beta  TRUE
## 4  8.2128315  8.5227315  8.0     h1 FALSE
## 5 11.6305653 11.9513433 12.0     h2 FALSE
## 
## $m8
##          low      high true params inint
## 1  3.9378608  4.557596  4.0    c.1  TRUE
## 2  5.7729105  6.447308  6.0    c.2  TRUE
## 3  0.4702384  0.533441  0.5   beta  TRUE
## 4  8.2315496  8.543833  8.0     h1 FALSE
## 5 11.6784857 11.993660 12.0     h2 FALSE
## 
## $m9
##          low       high true params inint
## 1  3.9731254  4.5753073  4.0    c.1  TRUE
## 2  5.6989545  6.3508008  6.0    c.2  TRUE
## 3  0.4805663  0.5419001  0.5   beta  TRUE
## 4  8.1834663  8.5009195  8.0     h1 FALSE
## 5 11.6774571 11.9984272 12.0     h2 FALSE
## 
## $m10
##         low      high true params inint
## 1 11.159384 12.494009   16    c.1 FALSE
## 2 18.014987 19.826380   24    c.2 FALSE
## 3  1.448803  1.602843    2   beta FALSE
## 4  7.672810  7.841917    8     h1 FALSE
## 5 12.192592 12.356752   12     h2 FALSE
## 
## $m11
##         low      high true params inint
## 1 11.185234 12.402412   16    c.1 FALSE
## 2 17.852498 19.520479   24    c.2 FALSE
## 3  1.445014  1.587621    2   beta FALSE
## 4  7.656411  7.817452    8     h1 FALSE
## 5 12.170768 12.333943   12     h2 FALSE
## 
## $m12
##         low      high true params inint
## 1 10.848001 12.049589   16    c.1 FALSE
## 2 17.692949 19.370861   24    c.2 FALSE
## 3  1.421928  1.564056    2   beta FALSE
## 4  7.664332  7.829884    8     h1 FALSE
## 5 12.197229 12.360644   12     h2 FALSE
## 
## $m13
##         low      high true params inint
## 1 12.524896 13.899434   16    c.1 FALSE
## 2 19.839749 21.731258   24    c.2 FALSE
## 3  1.604186  1.765372    2   beta FALSE
## 4  7.760157  7.913378    8     h1 FALSE
## 5 12.146768 12.297673   12     h2 FALSE
## 
## $m14
##         low      high true params inint
## 1 12.320277 13.634744   16    c.1 FALSE
## 2 19.447600 21.265906   24    c.2 FALSE
## 3  1.592442  1.745851    2   beta FALSE
## 4  7.735040  7.883800    8     h1 FALSE
## 5 12.148636 12.299303   12     h2 FALSE
## 
## $m15
##         low      high true params inint
## 1 12.213129 13.531782   16    c.1 FALSE
## 2 19.267255 21.137336   24    c.2 FALSE
## 3  1.575907  1.733495    2   beta FALSE
## 4  7.729892  7.877785    8     h1 FALSE
## 5 12.143976 12.292382   12     h2 FALSE
## 
## $m16
##         low      high true params inint
## 1 15.299240 16.740585   16    c.1  TRUE
## 2 23.408511 25.483513   24    c.2  TRUE
## 3  1.930559  2.106695    2   beta  TRUE
## 4  7.854474  7.987955    8     h1 FALSE
## 5 12.089570 12.224809   12     h2 FALSE
## 
## $m17
##         low      high true params inint
## 1 15.322033 16.795263   16    c.1  TRUE
## 2 23.497265 25.605068   24    c.2  TRUE
## 3  1.939864  2.115915    2   beta  TRUE
## 4  7.845753  7.980946    8     h1 FALSE
## 5 12.070039 12.206998   12     h2 FALSE
## 
## $m18
##         low      high true params inint
## 1 14.975978 16.394557   16    c.1  TRUE
## 2 23.018115 24.980390   24    c.2  TRUE
## 3  1.907683  2.074526    2   beta  TRUE
## 4  7.856500  7.991180    8     h1 FALSE
## 5 12.098570 12.233131   12     h2 FALSE

Gamma prior recaptures parameters relatively well here and there are no obvious model fitting problems. However, it pretty much always overestimates beta and cutpoints (while accurately identifying the inflection points). It’s not especially sensitive to the prior on \(\beta\) either. An \(\mathrm{exponential}\) prior with rates of 1,2,3, and 6 on \(\beta\) performed similarly.